//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"


// test for b field in the center and along axis of a circular current

#define RADIUS 1.0
#define CURRENT 5000
#define DIST1 0.1
#define DIST2 0.65
#define DIST3 3.65
#define DIST4 8.65
#define B_CENTER 0.0031415926535898
#define B_DIST1 0.0030950510166459
#define B_DIST2 0.0018517035538976
#define B_DIST3 0.0000579590991016
#define B_DIST4 0.0000047583098957210536
#define ALLOWED_FACTOR 1000.0 // allow difference by 1/1000 worth
#define Q 0

static void test_against (const char *where, const vect3d& bf, const vect3d& t_bf)
{
  prec_t diff = (bf - t_bf).length();
  if ( diff > t_bf.length() / ALLOWED_FACTOR )
    {
      printf ("ERROR_%s=fields differ too much\n", where);
    }
  else
    {
      printf ("OK_%s=fields match up\n", where);
    }
  char buf[128];
  printf ("BF_%s=%s\n", where, bf.sprint (buf));
  printf ("TBF_%s=%s\n", where, t_bf.sprint (buf));
  printf ("MARGIN_%s=%g (%g%%)\n", where, prec2double(diff), prec2double(100.0 * diff / t_bf.length()));
}

#if 0
#define DL 1e-5
#undef CURRENT
#define CURRENT 1
#undef Q
#define Q 1

struct mm : differator<4>
{
  mm (Statics& statics) : differator<4>(vect3d(DL, DL, DL)), statics(statics) { }
  Statics& statics;
  void getValues (const vect3d& pos, prec_t *vals)
  {
    vect3d bf, ef;
    statics.getFields (pos, bf, ef);
    vals[0] = bf.z;
    vals[1] = bf.x;
    vals[2] = ef.z;
    vals[3] = ef.x;
  }
};
#endif

int main (int argc, char *argv[])
{
  Statics statics;
  ElectroDynamics dyn(statics);
  StaticRing *ring = new StaticRing(CURRENT, vect3d(0, 0, 0), RADIUS, 0, vect3d(0.0, 0.0, 1.0));
  ring->setCharge (Q / (PREC_PI * RADIUS * 2));
  statics.addStaticElement (ring);

  vect3d bf;
#if 1
  dyn.getBField (vect3d(0, 0, 0), bf);
  test_against ("CENTER", bf, vect3d (0, 0, B_CENTER));
  dyn.getBField (vect3d(0, 0, DIST1), bf);
  test_against ("DIST1", bf, vect3d (0, 0, B_DIST1));
  dyn.getBField (vect3d(0, 0, DIST2), bf);
  test_against ("DIST2", bf, vect3d (0, 0, B_DIST2));
  dyn.getBField (vect3d(0, 0, DIST3), bf);
  test_against ("DIST3", bf, vect3d (0, 0, B_DIST3));
  dyn.getBField (vect3d(0, 0, DIST4), bf);
  test_against ("DIST4", bf, vect3d (0, 0, B_DIST4));
#endif

#if 0
  std::vector<vect3d> p;
  p.push_back (vect3d(0.1, 0, 3.1));
#if 1
  p.push_back (vect3d(0, 0, 0));
  p.push_back (vect3d(1e-5, 0, 0));
  p.push_back (vect3d(-1e-5, 0, 0));
  p.push_back (vect3d(0, 0, 0.8));
  p.push_back (vect3d(0, 0, -0.8));
  p.push_back (vect3d(0.8, 0, 0));
  p.push_back (vect3d(-0.8, 0, 0));
  p.push_back (vect3d(0.8, 0, 0.8));
  p.push_back (vect3d(0.8, 0, -0.8));
  p.push_back (vect3d(-0.8, 0, 0.8));
  p.push_back (vect3d(-0.8, 0, -0.8));
#endif

#if 0
  for ( size_t i = 0; i < p.size(); i++ )
    {
      vect3d pos = p[i];
      prec_t Ba, dBada, dBadr, d2Badadr, Br, dBrda, dBrdr, d2Brdadr;
      calcCoilBField (RADIUS, pos.z, pos.x, Ba, dBada, dBadr, d2Badadr, Br, dBrda, dBrdr, d2Brdadr);
      dyn.getBField (pos, bf);
      char b[128];
      printf ("a=%.3g, r=%.3g: a=(%.3g, %.3g, %.3g, %.3g) r=(%.3g, %.3g, %.3g, %.3g) a=%.3g, r=%.3g\n",
          prec2double(pos.z), prec2double (pos.x),
          prec2double(Ba), prec2double(dBada), prec2double(dBadr), prec2double(d2Badadr),
          prec2double(Br), prec2double(dBrda), prec2double(dBrdr), prec2double(d2Brdadr),
          prec2double(bf.z), prec2double(bf.x));

      mm d(statics);
      prec_t dfdx[4], dfdy[4], dfdz[4], d2fdxdy[4], d2fdxdz[4], d2fdydz[4], d3fdxdydz[4];
      d.getDifference (pos, dfdx, dfdy, dfdz, d2fdxdy, d2fdxdz, d2fdydz, d3fdxdydz);
      printf ("differator: a=(%.3g, %.3g, %.3g, %.3g) r=(%.3g, %.3g, %.3g, %.3g)\n",
          prec2double(bf.z), prec2double(dfdz[0]), prec2double(dfdx[0]), prec2double(d2fdxdz[0]),
          prec2double(bf.x), prec2double(dfdz[1]), prec2double(dfdx[1]), prec2double(d2fdxdz[1]));
    }
#endif

  for ( size_t i = 0; i < p.size(); i++ )
    {
      vect3d ef;
      vect3d pos = p[i];
      prec_t Ea, dEada, dEadr, d2Eadadr, Er, dErda, dErdr, d2Erdadr;
      calcCoilEField (RADIUS, pos.z, pos.x, Ea, dEada, dEadr, d2Eadadr, Er, dErda, dErdr, d2Erdadr);
      dyn.getEField (pos, ef);
      char b[128];
      printf ("a=%.3g, r=%.3g: a=(%.3g, %.3g, %.3g, %.3g) r=(%.3g, %.3g, %.3g, %.3g) a=%.3g, r=%.3g\n",
          prec2double(pos.z), prec2double (pos.x),
          prec2double(Ea), prec2double(dEada), prec2double(dEadr), prec2double(d2Eadadr),
          prec2double(Er), prec2double(dErda), prec2double(dErdr), prec2double(d2Erdadr),
          prec2double(ef.z), prec2double(ef.x));

      mm d(statics);
      prec_t dfdx[4], dfdy[4], dfdz[4], d2fdxdy[4], d2fdxdz[4], d2fdydz[4], d3fdxdydz[4];
      d.getDifference (pos, dfdx, dfdy, dfdz, d2fdxdy, d2fdxdz, d2fdydz, d3fdxdydz);
      printf ("differator: a=(%.3g, %.3g, %.3g, %.3g) r=(%.3g, %.3g, %.3g, %.3g)\n",
          prec2double(ef.z), prec2double(dfdz[2]), prec2double(dfdx[2]), prec2double(d2fdxdz[2]),
          prec2double(ef.x), prec2double(dfdz[3]), prec2double(dfdx[3]), prec2double(d2fdxdz[3]));
    }
#endif

  return 0;
}

